function out_t = rayleigh_t(fmax, fs, Ns)
N = 8;

while (N)
	if (N < 2*fmax*Ns/fs)
		N = 2*N;
	else
		break;
	end
end

N_inv = ceil(N*fs/(2*fmax));

f_interval = 2*fmax/N;
T_inv = 1/fs;

I_int = randn(1,N);
Q_int = randn(1,N);

I_inf = fft(I_int);
Q_inf = fft(Q_int);

sez(1) = 1.5/(pi*fmax);
for j = 2:N/2
	f(j) = (j-1)*f_interval;
	sez(j) = 1.5/(pi*fmax*sqrt(1-(f(j)/fmax)^2));
	sez(N-j+2) = sez(j);
end

%f = 0:f_interval:fmax-f_interval

%sez_temp = sez*(1/sqrt(1 - (f/fmax).^2));

k = 3;
p = polyfit(f(N/2-k:N/2),sez(N/2-k:N/2), k);
sez(N/2+1) = polyval(p, f(N/2)+f_interval);

I_outf = I_inf.*sqrt(sez);
Q_outf = Q_inf.*sqrt(sez);

I_temp = [I_outf(1:N/2) zeros(1, N_inv-N) I_outf(N/2+1:N)];
Q_temp = [Q_outf(1:N/2) zeros(1, N_inv-N) Q_outf(N/2+1:N)];

I_outt = ifft(I_temp);
Q_outt = ifft(Q_temp);

channel = I_outt + 1j* Q_outt;
channel = channel(1:Ns);

rms = sqrt(norm(channel, 'fro'));
channel_unit = channel/rms;

out_t = diag(channel_unit'*channel_unit);






